from EOS_collection import NRTL_MATRIX, antoine, gamma_cal, G_ex_MATRIX
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

#2024.08.05
#NRTL方程的拟合
#2024.12.11
#tau12,tau21与温度简单关联
#2025.1.22
#可以用于n个组分体系的拟合

# 常数
R = 8.314
atm = 101.325
Pa = atm
alpha = 0.3

# 基础参数
T = np.array([365.54, 364.88, 364.12, 363.19, 362.13, 361.21, 360.05, 359.04, 358.02, \
              357.53, 356.70, 355.56, 354.35, 353.91, 353.24, 352.32, 351.488, 350.76, \
              350.02, 348.61, 348.00, 347.21, 345.98, 344.33, 342.85, 340.56])
x_1 = np.array([0.0484, 0.0539, 0.0607, 0.0696, 0.0806, 0.0911, 0.1056, 0.1197, \
              0.1354, 0.1436, 0.1586, 0.1815, 0.2092, 0.2202, 0.2380, 0.2648, \
              0.2918, 0.3166, 0.3442, 0.4015, 0.4282, 0.4644, 0.5240, 0.6091, 0.6899, 0.8213])
y_1 = np.array([0.2704, 0.3009, 0.3213, 0.3492, 0.3866, 0.4126, 0.4492, 0.4762, \
              0.4885, 0.5047, 0.5323, 0.5677, 0.5969, 0.6031, 0.6237, 0.6519, \
              0.6681, 0.6923, 0.7197, 0.7431, 0.7443, 0.7613, 0.7860, 0.8231, 0.8633, 0.9221])
x_1_exp = x_1
y_1_exp = y_1
x = np.vstack((x_1, 1 - x_1))
y = np.vstack((y_1, 1 - y_1))
A = np.array([7.19736, 7.07404])
B = np.array([1574.99, 1657.46])
C = np.array([-34.29, -46.13])
dim = np.size(A)

# 控制输出的精度
np.set_printoptions(precision=3)

# 计算饱和蒸汽压
T1 = np.tile(T, (dim, 1)).T #在行上复制n次，并转置
Pvs = antoine(A, B, C, T1)

# 计算活度系数
gamma = gamma_cal(x.T, y.T, Pa, Pvs)
# print(f"gamma={gamma}")

def NRTL_fit(unknown, alpha):
    size = int(np.size(unknown) / 2) #拆分a，b两个参数
    a = unknown[ : size]
    b = unknown[size : ]
    gamma_cal = []
    for row_x, row_T in zip(x.T, T):
        tau_unknown = a + b / row_T
        tau = np.zeros((dim, dim))
        tau[np.where(np.eye(dim)==0)] = tau_unknown.flatten() #除了主对角线，其余依次赋值
        G = np.exp(-alpha * tau)
        gamma_temp = NRTL_MATRIX(row_x.reshape(-1, 1), tau, alpha) #将x变成列向量
        gamma_cal.append(gamma_temp)
    gamma_cal = np.vstack(gamma_cal) #按行堆叠
    f = np.sum( np.abs( gamma_cal - gamma ) ) ** 2
    return f

guess = np.ones( (dim * dim - dim) * 2 )  #tau的维度，一共n*n个，其中n个为0，参数a，b的数目为tau的2倍
result = minimize(NRTL_fit, guess, alpha, method = 'L-BFGS-B').x
a = result[ : dim]
b = result[dim :]

#计算相图
x_1 = np.linspace(0, 1, 20)
y_1 = np.ones_like(x_1)
x = np.vstack((x_1, 1 - x_1))
y = np.vstack((y_1, 1 - y_1))
T_guess = 300
Pa = atm

#求温度
def equilibrium(T, x):
    size = int(np.size(x))
    tau_known = a + b / T
    tau = np.zeros((size, size))
    tau[np.where(np.eye(size)==0)] = tau_known.flatten()
    G = np.exp(-alpha * tau)
    gamma = NRTL_MATRIX(x.reshape(-1, 1), tau, alpha)
    Pvs = antoine(A, B, C, T)
    Pa_cal = np.sum(Pvs * gamma * x)
    f = (Pa - Pa_cal) ** 2
    return np.sum(f) 

T2 = []
gamma = []
for row_x in x.T:
    size = int(np.size(row_x))
    T_temp = minimize(equilibrium, T_guess, args = row_x, method = 'L-BFGS-B').x
    tau_known = a + b / T_temp
    tau = np.zeros((size, size))
    tau[np.where(np.eye(size) == 0)] = tau_known.flatten()
    gamma_temp = NRTL_MATRIX(row_x.reshape(-1, 1), tau, alpha)
    T2.append(T_temp)
    gamma.append(gamma_temp)

T2 = np.vstack(T2)
gamma = np.vstack(gamma)
Pvs = antoine(A, B, C, T2)
y = Pvs * x.T * gamma / Pa
y_1 = y[:,0]

plt.rc('font',family='Times New Roman') #将全局的字体改为“Times New Roman”的形式
plt.rcParams['xtick.direction'] = 'in' #将坐标轴设置为朝内
plt.rcParams['ytick.direction'] = 'in' #将坐标轴设置为朝内
plt.plot(x_1,T2,label='NRTL')
plt.plot(y_1,T2,label='NRTL')
plt.plot(x_1_exp,T,'^',label='exp')
plt.plot(y_1_exp,T,'>',label='exp')
plt.xlabel('component 1 fraction')
plt.ylabel('T/K')
plt.legend()
plt.show()

